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Abstract 

A  laser-illuminated  imaging  system  operating  in  the  presence  of  atmospheric 
turbulence  will  encounter  several  sources  of  noise  and  diffraction  induced  errors.  As  the 
beam  propagates,  turbulence  induced  tilt  will  cause  the  beam  to  wander  off  axis  from  the 
target.  This  is  especially  troublesome  when  imaging  satellites,  since  most  turbulence  is 
closer  to  the  Earth’s  surface  and  greatly  affects  the  beam  in  the  early  stages  of 
propagation.  Additionally,  the  returning  beam  convolved  with  the  target  will  encounter 
turbulence  induced  tilt  that  appears  as  apparent  movement  of  the  target  between  image 
frames.  This  results  in  varying  beam  intensities  at  the  target  surface  between  imaging 
frames  that  can  affect  registration  algorithms  and  tracking.  In  this  research  effort,  an 
algorithm  using  expectation  maximization  and  least  squares  techniques  was  developed 
that  has  the  ability  to  separately  estimate  both  the  tilt  of  the  pulsed  laser  beam  and  the 
apparent  movement  of  the  object  between  incoherent  frames  and  produce  a  superior 
image  estimate  of  the  target  and  provide  tracking  information.  The  results  from  this 
algorithm  can  be  used  to  reduce  the  effects  of  beam  wander  and  increase  the  SNR  of 
post-processed  images. 
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Laser  Illuminated  Imaging: 

Multiframe  Beam  Tilt  Tracking  and 
Deconvolution  Algorithm 

I.  Introduction 

This  chapter  describes  the  problem  to  be  addressed  by  this  research,  the  goals  of 
this  project  and  previous  related  efforts.  Additionally,  the  assumptions  used  in  this 
research  on  the  system  and  its  data  are  examined  in  order  to  limit  the  scope  of  the 
problem.  Lastly,  an  outline  for  the  organization  of  this  thesis  is  given. 

1.1  Problem  Statement 

When  obtaining  high  resolution  images  from  a  laser  detection  and  ranging 
(LADAR)  system  there  are  several  factors  that  limit  the  systems  performance 
investigated  in  this  research,  such  as  diffraction  due  to  the  optics,  atmospheric  turbulence 
and  laser  beam  speckle.  These  factors  can  severely  distort  the  image  quality  and  reduce 
the  resolution  of  the  measured  data.  Due  to  operating  conditions  and  factors  such  as  cost, 
size  and  weight,  an  adaptive  optics  approach  may  not  be  feasible  for  all  situations.  This 
is  where  the  benefits  of  a  post-processing  algorithm  can  be  exploited  to  improve  the 
quality  of  LADAR  obtained  imagery.  The  proposed  algorithm  is  capable  of  providing 
estimates  of  the  target  image  with  distortions  such  as  speckle,  blurring  and  defocus 
mitigated  via  a  multiframe  processing  strategy. 

Atmospheric  turbulence  causes  random  time  delays  in  light  as  it  propagates 
through  the  atmosphere.  In  a  LADAR  application,  diffraction  due  to  the  atmospheric 
turbulence  results  in  tilt,  blur,  defocus  and  other  distortions  to  the  image.  Of  these 
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distortions,  tilt  accounts  for  87%  of  the  error  [14].  Using  a  Fourier  optics  approach  this 
time  delay  or  tilt  in  the  propagation  field  can  be  represented  as  a  spatial  shift  in  the  image 
field  [9].  Thus,  each  pulse  of  the  laser  beam  is  randomly  shifted  to  a  different  position  on 
the  target  and  the  returning  field  after  propagation  is  again  shifted  and  blurred.  This  shift 
in  the  image  field  is  what  a  registration  algorithm  corrects  for,  which  facilitates  averaging 
multiple  frames.  Additionally,  in  a  LADAR  system,  speckle  is  a  significant  source  of 
noise.  Speckle  is  caused  by  the  coherency  of  the  illuminating  laser  source  combined  with 
the  rough  surface  of  the  target  [10].  Each  frame  of  data  will  contain  independent 
intensity  fluctuations  that  appear  as  bright  and  dark  spots  as  a  result  of  laser  speckle. 
This  noise  further  complicates  the  registration  and  averaging  of  multiple  frames  of  an 
image  set. 

Current  deconvolution  and  registration  algorithms  are  successful  in  mitigating 
these  effects  however  they  do  not  take  into  account  the  effects  of  the  illuminating  beam 
shifting  around  due  to  turbulence  known  as  beam  wander.  Current  algorithms  assume  the 
beam  is  stationary  for  each  frame  of  data.  This  is  not  a  true  assumption  when 
propagating  a  beam  in  a  turbulent  atmosphere  with  a  beam  width  at  the  target  is  smaller 
than  the  field  of  view  (FOV)  of  the  receiver  optics. 

1.2  Research  Goals 

The  primary  goal  of  this  research  is  to  derive  and  prove  that  when  using  a 
multiframe  deconvolution  expectation  maximization  (EM)  algorithm  that  tracks  and 
estimates  the  beam  wander  in  each  frame  of  an  image  set,  the  global  shift  or  scene  shift 


2 


estimate  improves.  The  resulting  registered  image’s  spatial  resolution  will  improve  while 
also  providing  an  estimate  of  the  beam’s  position  at  each  frame. 

The  potential  increase  in  resolution  of  the  post-processed  image  and  the  tracking 
information  provided  has  both  commercial  and  defense  applications.  Such  applications 
include  laser  weapons  that  require  tracking  beam  wander  to  focus  a  high  energy  laser  on 
an  intended  target  and  remote  tracking  and  imaging  applications  that  involve  a  scene 
shifting  independently  from  a  shifting  beam. 

1.3  Assumptions 

For  this  research,  several  assumptions  of  the  situation  were  made  to  limit  the 
scope  of  the  project: 

•  The  shape  of  the  point  spread  function  (PSF)  due  to  the  receiver  optics  and 
atmospheric  turbulence  is  known  or  can  be  measured 

•  The  size  of  the  laser  beam  at  the  target  is  the  factor  limiting  the  FOV  of  the 
detected  image 

•  The  target  is  stationary  across  the  set  of  images  used  by  the  algorithm 

•  The  range  to  the  target  or  the  LADAR  beam  size  at  the  target  is  known 

•  The  mean  background  noise  is  known  or  can  be  measured 

•  The  beam  movement  is  small  enough  between  frames  so  that  the  portion  of 
the  target  being  illuminated  is  not  completely  different  than  the  other  image 
frames 

1.4  Thesis  Organization 

In  this  thesis,  Chapter  II  provides  the  background  and  theory  needed  to 
understand  the  concepts  in  this  research.  Additionally  Chapter  II  explains  the  effects  of 
beam  wander  on  current  image  registration  algorithms.  Chapter  III  describes  the 
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mathematical  derivation  of  the  proposed  algorithm,  along  with  the  simulation  and  the 
laboratory  set-up  developed  to  evaluate  the  algorithm.  Chapter  IV  details  the  results  and 
provides  an  analysis  from  testing  the  algorithm  using  the  experimental  set-up  and  the 
simulation  described  in  Chapter  III.  Lastly,  Chapter  V  gives  a  summary  of  the  research, 
provides  conclusions  on  the  thesis  and  offers  opportunities  for  future  work  to  expand  this 
effort. 
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II.  Background  &  Theory 


This  chapter  provides  an  overview  of  the  technical  background  material  necessary 
for  understanding  the  concepts  of  this  research.  A  brief  description  of  a  LADAR  system 
and  the  issues  that  affect  the  image  quality  and  the  ability  to  process  these  images  is 
provided.  A  brief  review  is  conducted  on  different  deconvolution  algorithms  that  have 
been  developed  and  their  limitations  to  the  addressed  problem.  The  effects  of  beam 
wander  on  the  most  common  image  registration  algorithm  is  illustrated  and  discussed. 
Additionally,  several  image  registration  techniques  that  are  quick  to  implement  and 
execute  are  discussed. 

2.1  LADAR  System  Model 

The  generic  LADAR  imaging  system  represented  in  this  research  interrogates  a 
target  using  a  coherent  pulsed  laser.  If  the  light  is  coherent,  it  is  assumed  that  the  phase 
of  the  laser  beam  is  constant  over  the  integration  time  of  the  camera.  Once  the  target  is 
illuminated  by  the  laser  beam,  the  beam  reflects  off  the  target  and  the  returning  pulse 
goes  through  the  LADAR  optics  systems  to  form  an  image.  This  basic  system 
description  is  illustrated  in  Figure  1.  As  discussed  later  in  this  chapter,  the  propagation 
path  between  the  LADAR  system  and  the  target  contains  atmospheric  turbulence  that 
adversely  affects  the  detected  image’s  quality. 
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Reflected  Beam 


Target 


LADAR 


Receiver  field  of  view 


Figure  1:  Basic  LADAR  system  concept 


Using  Fourier  optics,  a  basic  model  for  the  image  obtained  by  the  system,  i,  is 
described  in  Equation  1,  where  o  is  the  target  multiplied  by  the  illuminating  beam,  b, 
convolved  with,  h,  the  PSF.  In  this  equation,  the  detector  plane  coordinates  are 
represented  as  x  and  y  while  the  coordinate  system  at  the  target  is  in  the  z,  w  plane.  The 
two  reference  planes  are  considered  to  be  square  with  N  pixels  in  each  direction. 


N  N 


(l) 


z  w 


The  PSF  or  the  impulse  response  describes  the  response  of  the  system  to  a  point 
source  and  includes  the  effects  of  the  imaging  system  and  the  atmosphere.  This  research 
assumes  that  the  PSF  is  known  or  can  be  measured.  By  definition  the  optical  transfer 
function  (OTF)  abbreviated  as  H(fx,  fy)  is  simply  the  Fourier  transform  of  the  PSF.  The 
variables  fxand  fy  represent  the  spatial  frequencies  of  the  OTF  in  two  dimensions. 

There  are  several  scenarios  for  how  the  beam  interacts  with  the  target  as  shown  in 
Figure  2.  These  are  dependent  on  the  optics  of  the  FADAR  system  and  the  size  of  the 
beam  at  the  target.  It  is  possible  that  the  target  is  flood  illuminated,  meaning  the  beam  is 
larger  than  the  target  and  the  FOV  is  larger  than  the  target,  thus  the  entire  scene  is 
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illuminated  and  imaged.  Another  possibility  is  that  the  receiver  optics  FOV  is  limiting 
the  detected  image  thus  the  entire  detected  image  is  illuminated  but  only  a  portion  of  the 
scene  is  measured  by  the  optics.  The  third  possibility,  and  the  one  used  in  this  research, 
is  that  the  beam  size  at  the  target  is  both  smaller  than  the  actual  target  and  the  FOV  of  the 
receiver  optics.  Several  important  factors  limit  the  performance  of  the  LADAR  system 
during  this  process  including  atmospheric  turbulence,  noise  and  coherent  laser  speckle. 
Each  of  these  factors  is  discussed  in  detail  below. 
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(a) 


(b) 


Figure  2:  Detected  image  limiting  factors:  (a)  Flood  illuminated  limited,  (b)  Detector  FOV  limited. 

(c)  Beam  width  limited. 


2.2  Noise 

The  majority  of  noise  in  a  LADAR  system  can  be  attributed  to  one  of  the 
following  sources:  photon  counting  noise,  laser  speckle  and  background  noise  [13].  Each 
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of  these  noise  sources  additively  contribute  to  lowering  the  signal-to-noise  ratio  (SNR)  of 
the  detected  image  and  thus  affecting  the  performance  of  any  post-processing  algorithms. 

2.2.1  Photon  Counting  Noise 

Photon  noise  or  shot  noise  is  a  result  of  the  charge  coupled  device  (CCD)  array  in 
the  LADAR  system’s  camera  being  a  discrete  device  that  counts  the  number  of  photons 
that  arrive  at  each  element  of  the  array  [1],  The  distribution  of  the  photons  is  modeled  as 
a  Poisson  process  where  photons  arrive  at  random  intervals  with  a  mean  number  of 
photons,  X,  arriving  over  a  set  time  interval.  The  mean  is  also  a  random  process  that  will 
differ  from  frame  to  frame.  The  probability  of  r  photons  being  counted  at  a  pixel  spot  is 
given  by  the  Poisson  probability  mass  function  (PMF)  shown  in  Equation  2.  Due  to 
photon  noise,  the  same  target  imaged  at  two  different  times  will  not  have  the  same 
intensity  at  each  pixel  spot  creating  an  uncertainty  in  the  pixel  value  between  two 
different  frames  of  data. 


2.2.2  Speckle  Noise 

Speckle  noise  caused  by  the  coherent  nature  of  the  laser  in  the  LADAR  system  is 
usually  the  largest  source  of  noise  [13].  Laser  speckle  is  the  result  of  imaging  a  surface, 
rough  on  the  scale  of  an  optical  wavelength,  using  the  coherent  light  from  a  laser. 
Considering  that  the  optical  wavelength  in  a  LADAR  system  is  around  1pm,  the  variation 
in  the  surface  roughness  can  be  as  small  as  1pm  and  cause  significant  phase  change  in  the 
measured  light  at  the  detector  from  each  point  on  the  target.  As  illustrated  in  Ligure  3, 
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when  this  incoming  light  field  is  randomly  delayed  due  to  the  rough  surface,  their  sums 
can  interfere  with  one  another  causing  bright  and  dark  spots  [5].  Dark  spots  would  occur 
due  to  deconstructive  interference  when  the  phases  of  the  fields  cancel  each  other  out. 
Constructive  interference  resulting  in  bright  spots  would  occur  when  the  fields  arrive  in 
phase  and  increase  the  intensity  level.  According  to  work  by  Dr.  Goodman  [10],  speckle 
is  a  double  stochastic  process  in  that  the  intensity  fluctuations  follow  a  Gamma 
distribution  with  a  certain  number  of  degrees  of  freedom  related  to  the  coherency  of  the 
illuminating  light.  This  is  combined  with  Poisson  photon  counting,  resulting  in  a 
negative  binomial  distribution. 


Figure  3:  Illustration  of  phase  coherency  and  speckle  noise  due  to  a  rough  target. 

The  variance  of  the  laser  speckle,  shown  in  Equation  3,  represents  that  of  a 
negative  binomial  distribution  [13]  and  is  dependent  on  the  coherency  of  the  light  and  the 
expected  number  of  photons  received.  An  example  of  speckle  phenomenon  is  shown  in 
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Figure  4,  the  image  has  a  random  intensity  modulation  that  severely  distorts  the  image 
quality  when  compared  to  the  image  obtained  from  an  incoherent  light  source. 

^speckle  =  ^signal]  (l  +  (3) 

Where, 

•  E[AiSj5nai]is  the  expected  number  of  photons 

•  /J  is  the  coherency  factor  of  the  light,  1  =  fully  coherent,  oo  =  fully  incoherent 


(a)  (b) 

Figure  4:  Speckle  noise  illustration:  (a)  Incoherent  illuminated  image,  (b)  Coherent  illuminated 

image. 


Speckle  noise  can  be  mitigated  by  using  a  time  average  of  properly  registered 
images  [2],  Depending  on  the  coherence  of  the  light,  the  speckle  pattern  introduced  to 
each  image  can  be  especially  troublesome  when  registering  multiple  frames.  If  image 
quality  is  poor,  an  algorithm  might  not  properly  register  each  image  frame  and  will  blur 
the  averaged  result.  However,  a  deconvolution  algorithm  used  to  remove  the  effects  of 


11 


the  atmosphere  should  be  able  to  improve  the  quality  because  the  blurring  effect  from 
improperly  registered  frames  is  similar  to  the  blurring  effect  caused  by  the  atmosphere. 

2.2.3  Background  Noise 

Background  noise  is  the  result  of  any  light  or  signal  aside  from  the  illuminating 
beam  that  is  measured  by  the  detector  [13].  There  are  various  sources  of  background 
noise  and  many  are  dependent  on  the  situation,  they  include  but  are  not  limited  to  the  sun, 
starlight,  city  lights  and  the  laser  light  reflecting  off  other  surfaces  and  bouncing  back  to 
the  detector.  In  some  situations,  it  is  possible  to  estimate  the  amount  of  background  by 
taking  images  with  the  illuminating  laser  turned  off  to  get  a  mean  background  light  level. 

2.3  Atmospheric  Turbulence 

Turbulence  in  Earth’s  atmosphere  is  caused  by  random  variations  in  temperature 
and  air  motion  that  changes  the  refractive  index  of  the  air  [12].  As  optical  waves 
propagate  in  Earth’s  atmosphere,  the  wave  is  distorted  by  the  changes  in  the  refractive 
index  of  the  air  it  is  traveling  through  causing  phase  shifts  in  the  propagated  wavefront, 
this  is  illustrated  in  Figure  5. 
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Figure  5:  Effects  of  turbulence  on  a  propagated  wavefront  [11]. 


Using  Fourier  optics  and  the  shift  theorem  shown  in  Equation  4,  a  phase  shift 
represented  by  the  exponential  term  in  the  temporal  domain  introduces  a  spatial  offset  or 
translational  shift  (a,  b )  in  each  dimension  of  the  spatial  domain  [9].  The  terms  in  the 
equation,  g  and  G  represent  a  Fourier  transform  pair  with  G  being  the  Fourier  transform 
of  g. 

3{  g(x  —  a, y  —  b)}  =  G(fx,fy)e^2n^a+^  (4) 


The  total  variance  in  the  phase  of  the  wavefront  field  as  a  result  of  turbulence  is 
described  by  Equation  5,  where  D  is  the  diameter  of  the  receiver  aperture  and  ro  is  the 
coherence  diameter  or  Fried’s  seeing  parameter  [10].  The  coherence  diameter  or  seeing 
parameter  is  used  to  describe  the  optical  quality  of  the  atmosphere  and  is  typically  around 
5-10  cm  for  average  viewing  sites  and  up  to  20  cm  for  the  best  viewing  sites. 
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(5) 


The  phase  variance  due  to  tilt  [14],  in  one  axis  is  given  by  Equation  6,  this 
variance  is  doubled  when  looking  at  both  axis.  Nearly  87%  of  the  total  phase  variance  is 
a  result  of  tilt  with  image  distorting  effects  such  as  blurring  and  defocus  making  up  the 
rest.  Tilt  can  be  mitigated  by  an  accurate  image  registration  algorithm.  With  tilt 
removed,  the  residual  phase  variance,  given  in  Equation  7,  represents  the  higher  order 
image  distortions  such  as  defocus  that  are  corrected  for  using  a  deconvolution  algorithm. 
Shown  here  is  the  phase  variance  in  the  one  dimensional  0  direction.  These  equations 
would  be  identical  for  each  of  the  translation  shifts  axes. 


2 


0.448 


(6) 


O  tilt  removed 


0.134 


(7) 


The  atmospheric  turbulence  interacts  with  the  beam  causing  random  phase  delays 
and  results  in  a  beam  that  has  been  shifted  off  axis  from  the  intended  target.  The 
turbulence  will  also  introduce  distortions  in  the  beam  intensity  at  the  target  known  as 
beam  breathing  and  beam  scintillation  [8],  however  these  effects  were  not  studied  in  this 
research.  The  returning  field  is  affected  in  the  same  way  with  atmospheric  turbulence 
resulting  in  tilt,  blur  and  other  higher  order  distortions  on  the  returned  image.  Using  the 
Fourier  shift  theorem  previously  shown  in  Equation  4,  the  phase  delay  or  wavefront  tilt  in 
the  reflected  image  translates  to  a  global  spatial  shift  in  the  detected  image.  Taking  into 
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account  the  beam  shift  and  global  shift  introduced  by  the  turbulent  atmosphere,  an  update 
to  Equation  1,  which  was  an  expression  for  the  image  obtained  from  a  LADAR  system,  is 
shown  in  Equation  8.  In  this  equation,  ik,  is  the  kth  measured  data  frame  in  a  given  set. 
The  global  shifts  for  each  kth  frame  of  data  are  represented  in  the  PSF  as  translational 
shifts  ak  and  [ik  in  the  detector  plane  coordinate  system  x  and  y.  The  beam  shifts  are 
represented  as  yk  and  ek  in  the  target  plane  coordinate  system  z  and  w. 


ikix.y)  =  ^  ^  o(z,w)b(z  -  yk,w  -  £k)h(x  -  z  -  ak,y  -  w  -  /?fc) 


(8) 


Equation  8  implies  that  each  frame  of  data  obtained  from  an  image  set  with  a 
stationary  target  will  contain  both  a  beam  that  has  shifted  thus  illuminating  a  different 
location  on  the  target.  Additionally,  the  target  has  an  apparent  movement  from  the 
previous  frame  due  to  the  global  shift  introduced  in  that  frame.  This  is  illustrated  in 
Figure  6  with  simulated  sequential  data  frames  that  show  the  effects  of  turbulence 
induced  independent  shifting  on  the  illuminating  beam  and  the  global  scene. 
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(a)  (b) 

Figure  6:  Effects  of  turbulence  induced  tilt  on  beam  and  scene:  (a)  Frame  1.  (b)  Frame  2. 


The  amount  of  tilt  in  each  frame  is  not  necessarily  statistically  uncorrelated  and 
independent  from  the  previous  frame  or  time  instance.  There  is  a  degree  of  correlation  in 
the  amount  of  tilt  that  needs  to  be  considered  in  order  to  accurately  model  the  temporal 
characteristics  of  atmospheric  induced  tilt.  The  tilt  correlation  function  with  details 
found  in  [13]  is  dependent  on  the  characteristics  of  the  LADAR  system  such  as  the  time 
between  pulses  and  the  size  of  the  aperture,  as  well  as  the  degree  of  turbulence  and  the 
wind  speeds  at  the  imaging  site.  As  the  wind  speed  increases,  the  correlation  in  the 
atmospheric  turbulence  phase  screen  decreases  to  a  point  where  there  is  zero  tilt 
correlation  at  each  frame  of  data. 
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2.4  Image  Registration 

Image  registration  is  the  process  that  compares  and  aligns  multiple  images  by 
estimating  the  spatial  relationship  between  them.  Typically  these  relationships  could  be 
described  by  simple  translation  (horizontal  and  vertical  shifts),  rotational,  or  scaling 
differences.  Only  translation  shifts  were  evaluated  in  this  research.  Proper  registration  is 
vital  to  obtaining  averaged  images  that  have  a  sufficient  SNR  for  further  post-processing 
the  data  and  making  it  usable  for  many  applications. 

There  are  many  image  registration  techniques,  several  of  which  are  identified  and 
compared  in  [2].  Some  common  methods  include  cross-correlation,  directional  searching 
and  block  matching.  This  research  utilizes  a  version  of  the  directional  search  method  in 
which  the  sum  of  squared  errors  (SSE)  cost  function  is  iteratively  evaluated  in  four 
translational  shift  step  directions  (up,  down,  left  and  right  or  respectively  positive  y, 
negative  y,  negative  x  and  positive  x).  The  search  continues  until  the  SSE  cost  function  is 
minimized.  This  approach  is  susceptible  to  a  local  minimum  value,  but  in  this  application 
the  search  area  is  minimized  by  the  limitation  on  the  variance  of  the  beam  shift  and 
global  shift  as  previously  discussed. 

The  benefits  of  proper  image  registration  are  illustrated  in  Figure  7.  When 
multiple  frames  of  speckled  data  similar  to  that  shown  in  Figure  4  are  averaged  but  not 
registered  correctly,  the  speckle  noise  is  reduced;  however  the  image  is  significantly 
blurred  due  to  the  motion  between  frames  not  being  corrected.  If  the  images  are  perfectly 
registered  as  shown  in  Figure  8,  averaging  speckled  frames  of  data  removes  the  speckle 
noise  and  results  in  a  much  higher  resolution  image.  It  is  evident  that  proper  registration 
is  of  significant  importance  when  dealing  with  coherent  based  imaging  systems. 


17 


Figure  7:  Example  of  improper  image  registration  when  averaging  10  speckled  frames 


Figure  8:  Example  of  perfect  image  registration  when  averaging  10  speckled  frames 


Missing  from  these  algorithms  is  the  ability  to  register  multiple  frames  of  data 
while  tracking  and  incorporating  a  beam  that  has  shifted  positions  at  each  individual 
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frame.  To  illustrate  this,  a  cross-correlation  algorithm  was  used  to  register  multiple 
frames  of  data  [2],  In  one  scenario,  the  illuminating  beam  was  stationary  as  is  the 
assumed  case  in  all  the  algorithms  previously  mentioned.  In  another  scenario,  beam 
wander  was  introduced  to  each  frame.  The  results  in  Figure  9  and  Figure  10  show  that 
the  pixel  error  in  each  dimension  ( x  and  y)  was  considerably  greater  when  beam  wander 
was  present  compared  to  a  stationary  beam.  The  difference  in  the  error  between  these 
two  figures  is  likely  a  result  of  the  properties  of  the  simulated  target  and  not  a  result  of 
the  algorithm.  These  results  illustrate  the  purpose  of  this  research.  Preliminary 
simulations  show  that  beam  wander  causes  significant  registration  errors  in  common 
image  registration  techniques.  If  beam  wander  can  be  estimated  in  each  frame  and 
corrected  for,  its  effect  on  registration  can  possibly  be  decreased. 


Figure  9:  Cross  correlation  registration  error  with  and  without  beam  wander  in  the  x  direction 
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Figure  10:  Cross  correlation  registration  error  with  and  without  beam  wander  in  the  y  direction 


2.5  Deconvolution  Algorithms 

A  deconvolution  algorithm,  in  which  diffraction  and  distorting  effects  of  the 
atmosphere  and  optical  system  are  removed,  is  equally  important  to  obtaining  high 
resolution  and  quality  post-processed  images.  Deconvolution,  unlike  the  more  difficult 
blind  deconvolution,  assumes  these  affects  are  known  or  can  be  measured  through 
knowledge  of  the  PSF.  Several  of  the  most  common  algorithms  found  in  image 
processing  include  the  Ayers-Dainty  blind  deconvolution  technique  [7]  and  the  multi¬ 
frame  blind  deconvolution  (MFBD)  algorithm  [4].  The  MFBD  algorithm  is  an  iterative 
EM  approach  to  computing  the  maximum  likelihood  estimate  of  the  unknown 
parameters.  A  benefit  when  working  with  EM  algorithms  is  that  their  convergence  is 
assured  since  the  algorithm  is  guaranteed  to  increase  the  likelihood  function  at  each 
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iteration.  Again,  missing  from  both  these  algorithms  is  the  ability  to  track  a  wandering 
beam  in  each  frame  of  data  and  then  use  that  information  to  improve  the  deconvolution 
capability  of  the  algorithm. 
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III.  Methodology  &  Testing 


The  proposed  approach  to  reduce  the  effects  of  the  atmospheric  turbulence  and 
speckle  noise  of  imagery  obtained  from  a  LADAR  system  is  to  develop  an  EM  algorithm 
that  estimates  the  global  shift  and  the  beam  shift  in  each  frame  independently.  This 
chapter  describes  the  mathematics  in  developing  this  algorithm  and  the  implementation 
issues  with  the  EM  solution.  An  alternate  solution  is  proposed  based  on  proven 
algorithms  that  can  provide  superior  performance.  An  overview  of  the  simulated  and 
measured  data  is  given  along  with  criteria  that  will  be  used  to  test  if  the  proposed 
algorithm  results  in  greater  performance. 

Throughout  the  derivations  in  this  section,  all  equations  are  written  using  a  one 
dimensional  coordinate  system.  This  compresses  the  lengthy  equations  that  can  be  easily 
generalized  into  two  dimensions.  The  complete  two  dimensional  equations  are  given  in 
the  final  step  of  the  derivation. 

3.1  Expectation  Maximization 

The  EM  algorithm  is  an  iterative  approach  to  computing  the  maximum-likelihood 
estimate  of  the  unknown  parameters  in  a  given  data  set.  Similar  to  many  deconvolution 
algorithms,  the  EM  deconvolution  algorithm  proposed  in  this  research  is  derived  using 
Poisson  statistics.  Poisson  statistics  have  the  following  properties  that  make  working 
with  them  mathematically  simple  [6]: 

•  The  mean  of  a  Poisson  process  is  equal  to  its  variance 

•  The  sum  of  multiple  independent  Poisson  distributions  is  another  Poisson 
distribution  with  its  mean  equal  to  the  sum  of  the  means 
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Following  the  steps  outlined  by  Dempster,  Liard  and  Rubin  [3],  an  EM  algorithm 
is  proposed  to  iteratively  estimate  the  unknown  parameters  that  maximize  the  expected 
log-likelihood  function. 

3.1.1  Statistical  Model  for  the  Incomplete  Data 

Defining  a  statistical  model  for  the  incomplete  data,  which  is  the  observed  or 
measured  data,  is  the  first  step  in  formulating  an  EM  deconvolution  algorithm.  This 
model  is  defined  in  Equation  9.  It  is  known  that  the  incomplete  data,  d,  is  an  image  array 
of  independent  Poisson  random  variables  containing  the  true  target  image,  o,  multiplied 
by  the  beam,  b.  with  an  unknown  translational  shift,  y,  and  convolved  with  a  PSF,  h.  that 
contains  a  global  translational  shift,  a.  The  measured  background  radiation  is  represented 
as  B  in  Equation  9.  Notice  that  each  frame,  k,  contains  an  independent  shift  in  both  the 
beam  and  global  scene.  The  coordinate  system  used  in  the  derivation  of  this  algorithm  is 
shown  in  Figure  11. 

N 

E [40)]  =  ^  o(z)b(z  -  Yk)h(x  -  z-ak)  +  B{x )  (9) 
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Figure  11:  Defined  coordinate  system  in  the  detector  and  target  plane 


3.1.2  Define  the  Complete  Data 

The  complete  data,  d  ,  which  is  the  diffraction  and  noise  removed  desired  data  has 
to  be  statistically  related  to  the  incomplete  data  previously  defined  in  Equation  9.  The 
relationship  chosen  between  the  two  is  shown  below  in  Equation  10.  Unlike  the 
incomplete  data,  dk,  this  quantity  is  not  measureable  and  is  estimated  using  the  EM 
algorithm. 

N 

dk(x)  =  'YjTk(jx,z)  +  B(x)  (10) 

Z 

3.1.3  Statistical  Model  for  the  Complete  Data 

A  statistical  model  for  the  complete  data  is  defined  so  that  the  relationship 
expressed  in  Equation  10  produces  the  correct  statistical  model  for  the  incomplete  data. 
If  the  complete  data  is  considered  to  be  set  of  Poisson  random  variables,  the  incomplete 
data  can  be  related  to  the  complete  data  through  Equation  11.  Choosing  the  complete 
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data  to  be  Poisson  is  acceptable  because  the  sum  of  Poisson  random  variables  is  also  a 


Poisson  random  variable  [10]. 

E[dfc(x,z)]  =  o(z)b{z  —  yk)h{x  —  z  —  ak) 


(11) 


The  statistical  model  for  the  complete  data  is  verified  through  the  mathematical 
equivalence  expressions  in  Equation  12.  Starting  with  Equation  10,  the  expected  value  is 
taken  and  the  relationship  in  Equation  1 1  is  substituted  in.  This  result  is  equivalent  to  the 
relationship  defined  in  Equation  9. 


c(*)  =  I  dk(x,z )  +  B  {pc') 


E[4W1  =  E 


^  cLk{x,z )  +  B{x)  —  ^  E[dfc(x,z)]  +  B{x) 


(12) 


E[dk(x)]  =  ^  o{z)b(z  —  yk)h(x  —  z  —  ak)  +  B{x) 


3.1.4  Formulate  the  Complete  Data  Log-Likelihood 

Using  the  Poisson  PMF  previously  shown  in  Equation  2  and  applying  the  model 
for  the  complete  data  defined  in  Equation  11,  results  in  the  complete  data  likelihood 
expression: 


P  [dk(x,  z)] 


o{z)b{z  —  yk)h{x  —  z  —  ak)dk(jc,z^e  °(z)6(z  Y0h(x  z  “t) 
dk{x,  z)\ 


(13) 


The  next  step  is  to  expand  Equation  13  by  solving  for  all  pixel  points  ( x ,  z)  across 
the  image  set  containing  K  frames.  This  produces  the  joint  probability  for  the  complete 
data  likelihood  shown  in  Equation  14.  Due  to  the  independence  of  each  pixel  and  frame, 
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the  PMFs  for  each  pixel  point  and  frame  of  data  are  simply  the  product  of  the  individual 


PMFs  [6]. 


K  N  N 

P[dk(x,  z)  V  x,z ]  =nnn 

k  x  z 


o(z)6(z  -  YkMx  -Z-  akyk(x,z)e-o(z)Hz-Yk)h(x-z-ak) 
dk(x,z)\ 


(14) 


Following  the  EM  derivation  steps,  the  natural  log  of  Equation  14  is  taken  to  get 
the  log-likelihood  function,  L.  shown  in  Equation  15.  The  log-likelihood  function  is 
defined  by  the  unknown  parameters  in  the  complete  data  (o,  ak,  yk).  Greatly  simplifying 
this  equation  is  that  the  natural  log  of  the  product  operator  on  multiple  expressions  is  a 
much  simpler  summation  operator  on  each  expression  individually. 


K  N  N 


L  (o,ak,Yk)  =  III  cLk(x,  z)  In (o(z)fe(z  —  YkV l(x  ~  z  —  afc)) 

k  x  z 

—  o{z)b(z  —  Yu)h(x  —  z  —  ak )  —  ln(dfc  (x,z)l) 


(15) 


3.1.5  Expectation  of  the  Complete  Data  Log-Likelihood 

The  expectation  step  of  the  EM  algorithm  takes  the  conditional  expectation  of  the 
complete  data  log-likelihood  derived  in  Equation  15  when  given  the  incomplete  data  and 
previous  or  old  estimates  of  the  unknown  parameters  (o,  «/(,  yf}.  The  conditional 
expectation  is  calculated  and  shown  below  in  Equation  16. 
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=  E 


E[(L(o,  a,  y)  |  dold ,  dk  (x),  afd ,  ykd  )] 

K  N  N 

dk(x,  z)  In (o(z)£>(z  —  yR)h(x  —  z  —  ak))  —  o(z)£>(z  —  yk)h(x  —  z  —  afc) 

k  x  z 


-  ln(4  (x,  z) !)  |  dold ,  dk  (x),  ,  y£M 

K  N  N 

=  XXX  E[^(x’z)ln(°(z)fc(z  -  Yk)h{x  -  z  -  aJc))|ooW,dJc(x),<w,KtW] 

k  x  z 

-E[0(z)b(z-KS)h(x-z-a,)|ooW,dk(x),akoW-ykW]-E[ln(3;(x,z)!)|8oW,dk(x),<M,KkoW] 


(16) 


Due  to  the  linearity  of  the  conditional  expectation  function,  each  term  in  the 
conditional  expectation  in  Equation  16  can  be  evaluated  independently. 


First  Term 

The  first  term  in  the  conditional  expectation  can  be  simplified  by  moving  all  the 
terms  that  are  not  dependent  on  the  conditional  parameters  outside  of  the  expectation 
function.  The  second  step  in  simplifying  this  term  is  to  recognize  that  this  conditional 
expectation  is  related  to  the  binomial  PMF,  this  derivation  is  given  in  Appendix  A.  The 
final  form  for  the  first  term  is  shown  in  Equation  17.  The  term,  ikld,  given  in  Equation 
18,  is  the  estimate  of  the  image  based  on  the  past  estimates  of  the  unknown  parameters. 


E  [dk(x,z)  In  (o(z)b(z  -  yk)h{x  -z-  ak ))  |  dold,  dk(x),  akld,ykld] 

=  E[4(x,z)|ooW,c4(x),  a°kld,Ykd]  In (o(z)b(z  -  yfc)h(x  -  z  -  afc)) 
dold{z)b(z  —  ykld)h(x  —  z  —  a£w)d(x) 


(17) 


ikold{x) 


Tn(o(z)fe(z  —  yky i(x  —  z  —  ak)) 


Where, 


N 

ikold(x)  =  £  dold{z)b{z  -  Y°kld)l i(x  -  Z  -  a°kld ) 


Z 


(18) 
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Second  Term 


Similarly,  the  second  term  in  the  conditional  expectation  can  be  simplified.  This 
term  is  not  a  function  of  any  of  the  unknown  parameters  and  thus  all  the  terms  can  be 
brought  outside  of  the  expectation  function.  This  term  cannot  be  dropped,  even  though  it 
is  not  dependent  on  the  conditional  parameters.  Its  value  is  not  constant  and  will  change 
at  each  iteration  as  the  estimate  of  o  is  updated.  Simplifying  leads  to  the  final  form 
shown  in  Equation  19. 

E[o(z)t>0  -  Yk)Kx  -z-  ak)\80ld,  dk(x),  a°kld,  yk‘ d]  =  o(z)b(z  -  yk)/i(x  -  z  -  ak)  (i9) 

Third  Term 

The  final  term  in  the  conditional  expectation,  containing  the  factorial  operation, 
can  thankfully  be  ignored.  The  term  is  not  conditional  on  old  estimates  of  the  unknown 
parameters  that  are  being  estimated  or  the  incomplete  data.  Thus  this  term  is  a  constant 
value  when  maximizing  the  function  with  respect  to  the  unknown  parameters  and  does 
not  need  to  be  evaluated  to  maximize  the  complete  data  log-likelihood  function. 

Total  Conditional  Expectation 

Combining  the  results  of  the  three  individual  pieces  of  the  conditional 
expectation,  leads  to  the  final  form  of  the  complete  data  log-likelihood  conditional 
expectation,  Q,  shown  in  Equation  20. 
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(20) 


Q  =  E [(L (o,  a,  y)  |  oold ,  dk  (x),  akld ,  ykld  )] 


■m 

k  x  z 
K  N  N 

III 


QOld  (z)b(z  _  y°ld  ~)h(x  ' 


1)dk(pc) 


i*oW0O 


gold  (z)b(z  _  y°ld  ~)h{x  ■ 


1)dk(x ) 


k  x  z 

-  o(z)b(z  -  yk)h(x  -z-ak) 


ik°ld(x ) 


In (o(z)h(z  -  yi:)h(x  -  z  -  ak))  -  o{z)b{z  -  yk)h(x  -  z  -  ak) 


(ln(o(z))  +  ln(fe(z  -  yk))  +  In  {h{x  -z-  ak ))) 


3.1.6  Maximization  of  the  Complete  Data  Log-Likelihood  Conditional  Expectation 

With  the  conditional  expectation  known  (Equation  20),  the  next  step  is  to 
maximize  it  with  respect  to  the  parameters  being  estimated.  This  process  is  completed 
separately  for  the  three  different  sets  of  unknown  parameters:  the  global  shift,  the  beam 
shift,  and  the  true  object.  In  each  instance,  the  terms  not  dependent  on  the  specific 
parameter  being  estimated  can  be  dropped  since  they  do  not  influence  maximizing  the 
function.  The  derivation  for  each  of  the  three  unknown  parameter  sets  is  given  below 
individually. 

Maximize  Global  Shift 

When  maximizing  the  conditional  expectation  from  Equation  20  with  respect  to 
the  global  shift  parameter  ( ak ),  the  terms  that  are  not  dependent  on  this  parameter  can  be 
dropped  as  shown  in  Equation  21.  Recognizing  that  the  summed  PSF  term  is  constant  for 
all  values  of  shifts  means  that  term  can  also  be  dropped.  Since  each  frame  contains  a 
different  shift,  each  frame  is  evaluated  independently  by  setting  k  —  k0. 
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Qa 


iV  iV 

III 

k=kn  x  z 


dold(z)b{z  —  ykld)h(x  -  z  -  akld}cLk(x) 


i*  (*) 


(in  (o(z)')  +  In  (b{z — y^} 


+  In  (h(. 


x  —  z  —  a 


fc))  )  -  o(z)b(z  —  yk)h(x  —  z — et^ 


(21) 


-II 


dold(z)b{z  -  ykld)h(x  -z-  a^ld)dko(x) 


iko°ldM 


In  (li(x-z-afco)) 


Thus,  the  final  form  of  the  maximum  likelihood  expression  for  the  global  shifts, 
Qa,  in  a  single  given  frame,  ko,  after  slight  rearranging  of  terms  is  shown  in  Equation  22. 


Qako  =  ^  d°ld(z)b(z  -  y£ld)  ^  ■  Mx  -  2  -  <d)ln  ( Kx  -  2  ~  «fc0)) 


‘k„oW00 


(22) 


Maximize  Beam  Shift 

When  maximizing  the  conditional  expectation  from  Equation  20  with  respect  to 
the  beam  shift  parameter  (yfc),  the  terms  that  are  not  dependent  on  this  parameter  can  be 
dropped  with  the  final  form  shown  in  Equation  23  having  been  slightly  rearranged. 
Again,  since  each  frame  contains  an  independent  shift,  each  frame  is  evaluated 
individually  by  setting  k  —  k0. 
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Lrkl  = 


.-III 
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Maximize  Target 

When  maximizing  the  conditional  expectation  from  Equation  20  with  respect  to 
the  target  (o),  the  process  is  slightly  different  than  the  previous  steps.  To  find  the 
maximum  likelihood  estimate,  the  derivative  of  the  conditional  expectation  with  respect 
to  a  single  pixel  point  z0  in  the  true  target,  o(z0),  is  calculated.  The  results  are  shown  in 
Equation  24. 
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Solving  for  the  derivative  of  the  first  term  of  Equation  24  gives  the  results  shown 
in  Equation  25. 
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Additionally,  the  derivative  of  the  second  term  of  Equation  24  is  shown  in 
Equation  26. 
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Combining  the  two  terms  in  Equation  25  and  Equation  26,  and  then  applying  the 
sifting  property,  results  in  the  final  form  of  the  conditional  expectation  that  will  be 
maximized  shown  in  Equation  27. 
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Equation  27  is  set  to  zero  to  find  when  the  maximum  value  will  occur.  This  result 
is  shown  in  Equation  28. 
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A  rearrangement  of  terms  in  Equation  28  and  simplifying  the  expression  by 


removing  the  summed  PSF  results  in  Equation  29. 
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Solving  for  the  true  target,  o(z0),  produces  the  final  update  equation  shown  in 
Equation  30.  At  each  iteration  of  the  algorithm  the  new  estimate  of  the  true  target  is 
identified  as  onew(z0).  This  solution  is  ideal  because  it  could  be  used  iteratively  and 
updates  the  target  at  each  iteration  based  only  on  past  estimates  of  the  unknown 
parameters  (akld,  kld,  o0,d).  Of  possible  concern  is  the  summation  of  all  the  shifted 
beams  in  the  denominator  of  Equation  30.  However,  mathematically  a  Gaussian 
approximation  to  long  term  beam  wander  as  a  result  of  atmospheric  turbulence  can  be 
taken  that  eliminates  the  requirement  to  sum  all  of  the  beams  individually  [12]. 


Lik  2-tx  ' 


ii.z0)b{z 


y°uld)h{.x 


vold 


)dk(x) 


7(z0)  = 


lw 


Ykb(z0  ~Yk) 


(30) 


Expectation  Maximization  Update  Equations 

Expanding  Equation  22,  Equation  23  and  Equation  30  to  their  two  dimensional 
version  results  in  the  EM  algorithm  solution  shown  in  Equation  31,  Equation  32  and 
Equation  33. 
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Examining  these  equations  shows  that  a  solution  to  solving  for  each  of  the 
parameters  individually  is  mathematically  possible  by  using  an  iterative  process  that 
would  first  update  the  object  estimate,  onew,  then  estimates  the  global  shift  parameters 
for  each  frame,  ak  and  /?fc,  and  then  estimates  the  beam  shifts,  Yk  and  sk.  However,  a 
tractable  solution  could  not  be  found  in  MATLAB  during  the  implementation  of  the  two 
maximum  likelihood  expressions  solving  for  the  beam  and  global  shifts.  The  likely  cause 
is  that  in  both  instances,  MATLAB  was  unable  to  evaluate  the  natural  log  of  the  slightly 
shifting  beam  or  PSF  with  enough  accuracy  to  correctly  estimate  the  shift.  The  non¬ 
linear  properties  of  the  natural  log  function  and  the  fact  that  the  beam  and  PSF  approach 
zero  on  the  tails  caused  the  change  between  the  beam  or  the  PSF  shifts  from  each  frame 
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to  be  smaller  than  the  machine’s  precision.  It  is  possible  that  the  algorithm  might  have 
converged  to  a  solution  if  given  more  time  or  a  computer  with  higher  precision. 

Although  the  pure  EM  approach  did  not  completely  work,  the  estimate  for  the 
target  update  worked  and  the  fact  that  an  EM  solution  exists  for  finding  the  shift  values 
based  only  on  past  estimates  suggest  that  solving  for  the  target  and  the  shifts  independent 
is  mathematically  possible.  To  attempt  to  solve  for  the  beam  and  global  shifts  a  different 
cost  function  was  studied. 

3.2  Two  Dimensional  SSE  Approach 

Along  with  the  target  update  from  Equation  30,  an  iterative  least  squares 
likelihood  cost  function  [2]  was  taken  to  solve  for  the  global  and  beam  shifts.  The  least 
squares  algorithm  steps  the  unknown  parameter  in  each  direction  (±a,  /?  or  +  y,  s) 
separately  and  calculates  the  error  at  each  location  using  Equation  34.  The  step  that 
results  in  the  least  error  is  the  direction  to  move  the  unknown  parameter  to. 

N  N 

SSE  =  Z  Z  K  O’  y)  -  <  O’  yj)  (34) 

x  y 

As  the  algorithm  iterates,  each  step  brings  the  estimate  closer  to  the  actual 
solution.  This  technique  is  accomplished  one  frame  at  a  time  for  first  the  global  shifts 
and  then  for  the  beam  shifts.  When  the  minimum  error  is  at  the  current  location,  the 
algorithm  stops  iterating.  Figure  12  shows  the  flow  of  this  technique  for  estimating  the 
global  shifts.  Estimating  the  beam  shift  is  not  shown  since  it  is  exactly  the  same  with  y 
and  s  substituted  in  for  a  and  p.  The  minimum  step  size  at  each  iteration,  SS,  is  defined 
in  the  algorithm  and  can  be  adjusted  based  on  computation  and  accuracy  requirements. 
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Figure  12:  Two  dimensional  SSE  algorithm  flow 


3.3  Testing  Methods 

The  algorithm  develop  in  this  research  was  tested  using  both  purely  simulated 
data  as  well  as  a  hybrid  measured  data  set.  These  data  sets  were  used  to  test  and  evaluate 
the  performance  of  the  proposed  algorithm  compared  to  current  registration  and 
deconvolution  algorithms. 


3.3.1  Simulated  Data 

The  simulated  images  were  generated  in  MATLAB  for  testing  the  proposed 
algorithm’s  ability  to  estimate  the  global  and  beam  shifts  while  deconvoluting  and 
registering  the  frames  to  obtain  an  estimate  of  the  target.  The  simulated  data  was  a  set  of 
30  image  frames.  The  target,  previously  shown  in  Figure  4,  is  a  512x512  pixel  array  of  a 
mobile  United  States  Air  Force  (USAF)  resolution  target  board.  The  illuminating  beam 
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was  simulated  to  have  a  two  dimensional  unit  Gaussian  intensity  profile  with  a  standard 


deviation  of  35  pixels  as  shown  in  Figure  13. 


Pixels  °  0  Pixels 

Figure  13:  Simulated  unit  Gaussian  illuminating  beam 


The  simulated  PSF  contained  a  generic  defocus  error  to  mimic  blurring  of  the 
target  due  to  atmospheric  turbulence.  To  simplify  the  problem,  but  still  capture  89%  of 
phase  aberrations  [14]  only  focus  and  tilt  errors  were  included  in  the  simulation. 

The  beam  and  global  shifts  in  each  frame  due  to  atmospheric  turbulence  were 
chosen  from  the  Gaussian  distribution  as  zero  mean  with  a  standard  deviation  of  3  pixels 
for  the  global  shift  and  a  7  pixels  standard  deviation  for  the  beam  using  the  Gaussian 
number  generator  function,  rand,  in  MATLAB.  This  allowed  enough  movement  in  both 
the  beam  and  scene  to  test  the  algorithm  without  entirely  changing  the  detected  image 
between  frames.  A  wind  velocity  of  10  m/s  across  the  aperture  of  the  camera  was 
assumed,  this  resulted  in  zero  tilt  correlation  between  frames  of  data  and  so  the  tilt  at 
each  frame  is  completely  independent  and  uncorrelated  with  the  previous  frame  [13]. 
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There  are  two  ways  that  speckle  could  be  simulated.  In  the  first  scenario,  a 
random  phase  screen  is  multiplied  with  the  product  of  the  beam  and  the  target.  This 
represents  the  phase  delay  expected  in  the  reflected  image  due  to  the  rough  surface  of  the 
target.  Then  using  Fourier  transforms  to  simulate  propagation,  the  phase  delay  resulted 
in  speckle  at  the  image  plane.  However,  an  alternate  approach  was  used  when  creating 
simulated  data  for  this  research.  First,  an  image  was  simulated  in  MATLAB  assuming 
completely  incoherent  light  was  being  used.  Then  speckle  noise  was  added  to  the 
predicted  image  at  the  detector  plane  using  the  MATLAB  icdf  function  to  add  negative 
binomial  noise  with  30  degrees  of  freedom  to  each  frame  to  represent  the  speckle 
expected  from  partially  coherent  light.  The  icdf  function  in  MATLAB  adds  the  desired 
random  distribution  to  the  input  variable  when  given  the  mean  value  and  a  certain 
number  of  degrees  of  freedom,  M ,  related  to  how  coherent  the  light  is  which  effects  the 
degree  of  intensity  fluctuations  [10].  This  process  simulates  speckled  data  that  contains 
the  same  statistics  as  would  be  expected  using  the  correct  physical  model  described  in  the 
first  scenario  since  the  expected  value  of  multiple  speckled  images  that  are  averaged  is 
the  image  expected  if  using  an  incoherent  light  source. 

The  background  noise  was  generated  using  the  Poisson  number  generator  with  a 
mean  of  20  photons  and  added  to  the  detected  image.  The  parameters  used  in  the 
simulated  data  are  summarized  in  Table  1. 
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Table  1  Simulation  setup  parameters 


Parameter 

Value 

Image  size 

512x512  pixels 

Beam  width  standard  deviation 

35  pixels 

Global  shift  standard  deviation 

3  pixels 

Beam  shift  standard  deviation 

7  pixels 

Aperture  diameter,  D 

2  mm 

Time  between  pulses 

.1  s 

Wind  velocity  across  aperture 

10  m/s 

Intensity  degrees  of  freedom,  M 

30 

Mean  background  noise,  B 

20  photons 

Detected  image,  max  photon  count 

250  photons 

Number  of  frames  in  data  set 

30  frames 

Shown  in  Figure  14  is  an  example  of  four  simulated  consecutive  frames  of  data. 
The  different  beam  and  global  shifts  are  apparent  between  frames  as  well  as  the 
distinctive  speckle  pattern  among  each  frame. 
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Frame  1 


Frame  2 


Figure  14:  Four  frames  of  simulated  data 
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3.3.2  Measured  Data 


It  was  unfeasible  to  collect  measured  data  that  contained  controlled  beam  shifting 
as  well  as  global  shifting  in  the  laboratory  environment  available  for  this  research.  To 
compensate,  a  hybrid  approach  was  taken  to  collect  measured  data.  An  imaging  system 
was  set  up  in  the  laboratory  with  the  computer  display  screen  at  the  focal  plane  of  a 
focusing  lens,  4.5  inches  in  front  of  the  camera  as  shown  in  Figure  15.  The  imaging 
system  captured  20  frames  of  the  same  simulated  data  frame  and  then  the  next  frame  of 
data  was  displayed  on  the  computer  screen.  The  20  frames  of  data  were  averaged  to 
obtain  a  higher  SNR  on  the  detected  image. 


2mm  Limiting  Aperture 


4.5"  Focusing  Lens 


97"  -  Propagation  Distance 


Computer 

Display 


Figure  15:  Laboratory  setup 


The  data  displayed  on  the  computer  screen  contained  the  target  multiplied  by  a 
shifted  beam  and  that  result  shifted  again  to  simulate  the  scene  shift.  This  was  generated 
using  MATLAB  in  a  similar  manner  as  was  described  in  the  simulated  data  with  several 
variations.  First,  the  target  multiplied  by  the  beam  was  not  convolved  with  a  PSF,  the 
convolution  with  the  PSF  occurred  naturally  in  the  setup  due  to  the  2mm  aperture. 
Second,  the  negative  binomial  distribution  of  speckle  was  simulated  by  adding  Gamma 
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distributed  noise  using  the  MATLAB  icdf  function.  Gamma  distributed  noise  was  used 
to  simulate  the  intensity  fluctuation  expected  due  to  the  rough  surface  target  [10].  The 
naturally  occurring  Poisson  process  (camera  photon  counting)  with  a  mean  that  has  a 
Gamma  distribution,  results  in  a  negative  binomial  distribution  expected  for  speckle 
noise.  Additionally,  background  noise  was  not  added  to  the  image  since  this  will  occur 
naturally.  Shown  in  Figure  16  are  the  first  four  frames  of  measured  data  after  20  frames 
of  identical  data  are  averaged  to  improve  the  SNR  of  the  frames  to  be  used  by  the 
algorithm. 
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Frame  1 


Frame  2 


Frame  3 


Frame  4 


Figure  16:  Four  frames  of  measured  data 

Collecting  measured  data  using  this  hybrid  approach  has  the  key  advantage  in  that 
the  truth  data  is  known  because  the  beam  and  global  shifts  were  defined  in  MATLAB. 
This  approach  can  still  be  called  measured  data  because  an  optical  system  is  used  to 
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capture  the  data,  the  speckle  noise  and  background  noise  is  truly  random  and  the  PSF  is 
convoluted  naturally  in  this  set-up.  However,  the  value  for  all  parameters  will  not  be 
known  exactly  as  is  the  case  in  a  simulated  environment.  The  parameters  for  the 
measured  data  set  up  are  summarized  in  Table  2. 


Table  2  Measured  data  parameters 


Parameter 

Value 

Width  of  camera  captured  image 

512  x  512  pixels 

Display  window  physical  size 

19.5  x  19.5  cm 

Display  size,  N 

512  pixels 

Beam  width  standard  deviation 

35  pixels 

Beam  shift  standard  deviation 

7  pixels 

Global  shift  standard  deviation 

3  pixels 

Coherence  parameter,  M 

30 

Number  of  identical  frames  averaged 

20 

Camera  integration  time 

0.1  sec 

Aperture  diameter,  D 

2  mm 

Focusing  lens  focal  length,  it 

4.5  in 

Distance  from  lens  to  display 

97  in 

Number  of  frames  in  data  set 

30  frames 

Pixel  pitch  on  camera 

16  um 

Background  radiation,  B  (measured) 

1036  photons 

Knowing  the  relationship  between  a  pixel  represented  on  the  computer  display 
and  a  pixel  captured  by  the  camera  system  is  an  important  piece  of  information  in  this 
hybrid  approach.  Following  principles  of  optics  systems  [9],  the  optics  magnification 
factor  can  be  calculated  using  the  property  of  similar  triangles  with  the  measurements 
shown  in  Figure  15  and  Table  2.  This  information  allows  a  shift  in  the  MATLAB 
environment  to  be  translated  into  a  shift  in  the  detected  image.  Shown  in  Equation  35, 
the  magnification  ratio,  Ar,  is  calculated  as  1.104.  Thus  a  shift  of  1  pixel  in  MATLAB 
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and  displayed  on  the  computer  screen  would  correspond  to  a  shift  of  1.104  pixels  in  the 


image  detected  by  the  camera. 


A,.= 


size  of  pixel  on  display 

- - - - - x  optics  magnification  = 

size  of  pixel  on  CCD 


(19.5  cm\ 
l  512  J 
16  um 


1.104 


(35) 


3.3.3  Comparison  Criteria 

Several  criteria  were  used  to  evaluate  the  performance  of  the  algorithm 
using  the  measured  and  simulated  data.  For  both  the  simulated  and  measured  data,  the 
error  of  the  estimated  shift  parameters,  E  ’,  at  each  frame  of  data  can  be  calculated  using 
Equation  36  when  the  true  shifts  are  known.  Additionally,  this  error  can  be  calculated  for 
both  scenarios  where  the  beam  is  being  tracked  in  the  algorithm  at  each  frame  and  when 
the  beam  is  considered  to  be  stationary.  To  be  successful,  the  proposed  algorithm  needs 
to  provide  a  better  estimate  of  the  global  shifts  when  the  beam  position  is  tracked.  This 
would  result  in  a  decrease  in  the  registration  error  and  thus  an  image  with  a  higher  SNR 
and  greater  resolution  when  multiple  frames  are  averaged. 

I’a  (.k)  r^frue(/c)  oc{kk) 

EpW  =  PtrueiX)  -  Kk) 

E^k)  =  Ytrue{k)-Y{k)  M(>! 

£'(fc)  =  £true(fe)  -  £(fc) 
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Using  the  simulated  data,  the  estimated  image,  d(z,  w,  k ),  at  each  iteration  k,  can 
be  compared  to  the  known  true  target,  o(z,  w),  to  calculate  the  root  mean  squared  error 
(RMSE)  at  each  iteration  using  Equation  37  where  N  is  the  number  of  pixels  in  the  image 
array.  The  RMSE  can  be  calculated  using  both  the  scenario  where  the  beam  is  tracked 
and  with  the  beam  tracking  off.  A  marked  improvement  in  RMSE  when  the  beam  is 
tracked  would  indicate  the  algorithm  provides  better  performance  than  standard  image 
registration  algorithms  that  do  not  track  beam  wander. 


RMSE{k )  = 


Zw=iEz=l  o(z,w,  k)  -  o(z,w) 


N2 


(37) 
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IV.  Results  and  Analysis 


This  chapter  presents  the  results  of  applying  the  algorithm  developed  in  Chapter 
III  to  both  the  simulated  and  measured  data  sets.  In  both  situations,  the  error  in  the 
estimated  shift  parameter  is  evaluated.  In  the  simulated  environment,  the  RMSE  of  the 
estimated  image  is  calculated  at  each  iteration.  These  metrics  are  analyzed  to  determine 
if  the  proposed  algorithm  provides  an  improved  performance  over  cross-correlation, 
which  is  the  most  commonly  used  image  registration  technique. 

4.1  Simulated  Data  Results 

Using  simulated  data  created  as  described  in  Chapter  III,  the  proposed  algorithm 
was  used  to  estimate  the  shifts  and  true  target  using  the  parameters  shown  in  Table  3. 
The  scale  or  step  size  for  estimating  the  shifts  was  set  to  a  quarter  of  a  pixel.  Changing 
the  step  size  allows  for  a  tradeoff  between  better  estimates  of  the  shifts  but  with  a 
significantly  longer  execution  time  as  the  step  size  is  decreased. 


Table  3  Algorithm  parameters 


Parameter 

Value 

Number  of  frames  of  data 

30 

Pixel  shift  scale 

!4  pixel 

Max  number  of  iterations 

50 

Using  simulated  data  with  the  algorithm  tracking  the  beam  produced  an  estimated 
image  of  the  target  shown  in  Figure  17  after  50  iterations.  The  exercise  was  repeated 
using  the  same  data  set  and  parameters  except  the  ability  to  track  the  beam  disabled.  This 
resulted  in  the  estimated  image  shown  in  Figure  18.  Visually  these  two  results  look 
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similar  with  the  exception  that  the  estimated  image  when  tracking  the  beam  is  slightly 
larger  than  the  estimate  obtained  when  the  beam  is  not  tracked.  This  is  due  to  the  beam 
being  estimated  at  each  frame  and  data  on  the  edges  of  each  frame  not  being  lost  when 
multiple  frames  are  averaged. 


Figure  17:  Estimated  target  using  simulated  data  -  beam  tracking  on 
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Figure  18:  Estimated  target  using  simulated  data  -  beam  tracking  off 


The  RMSE  of  the  estimated  image  after  each  iteration  was  calculated  using 
Equation  37  and  the  results  are  shown  in  Figure  19.  After  the  50th  iteration,  the  RMSE 
was  27.9  photons  with  beam  tracking  on  and  32.4  with  beam  tracking  off.  This 
represents  a  13.8%  improvement  in  RMSE  performance  when  the  beam  is  tracked  at  each 
frame. 
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Figure  19:  RMSE  -  Simulated  data  photon  error  at  each  iteration 


The  estimated  shifts  from  the  algorithm,  given  explicitly  in  Appendix  B  along 
with  the  true  shifts,  are  summarized  by  using  Equation  36  to  compute  the  error.  The 
calculated  error  in  the  estimated  shift  for  each  frame  of  data  when  compared  to  the  true 
shifts  is  shown  for  each  parameter  (a,  p,  y  and  s)  in  Figure  20,  Figure  21,  Figure  22  and 
Figure  23.  When  examining  the  beam  shift  error  in  Figure  22  and  Figure  23  there  is  only 
one  set  of  data  plotted  because  when  the  algorithm  is  not  estimating  the  beam  shift,  the 
error  of  the  estimated  beams  position  is  irrelevant. 
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Figure  20:  Error  in  shift  estimate  for  a  using  simulated  data 


Figure  21:  Error  in  shift  estimate  for  p  using  simulated  data 
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Analyzing  the  information  provided  in  the  previous  graphs  show  that  when  using 


simulated  data,  the  algorithm  provides  superior  results  when  the  beam  tracking  feature  is 
enabled.  The  difference  in  the  algorithm’s  ability  to  track  better  in  one  direction  versus 
the  other  is  likely  a  result  of  the  physical  properties  of  the  target.  One  concern  is  the 
outlier  in  frame  15  that  appears  in  Figure  21  and  Figure  22.  The  global  shift  estimate  for 
this  frame  is  wildly  off  from  all  other  estimates.  This  is  likely  a  result  of  the  combination 
of  a  large  beam  shift  and  the  fact  that  the  bar  charts  on  the  target  at  many  locations  look 
nearly  identical.  Due  to  the  beam  illuminating  a  different  bar  chart  on  the  target  in  that 
frame  caused  the  global  shift  to  be  off.  To  ensure  this  wasn’t  skewing  the  results,  the 
mean  of  the  error  in  the  estimated  shift  parameters  is  summarized  in  Table  4  with  the  15th 
frame  removed  from  the  mean  calculations.  When  using  the  proposed  algorithm  to 
estimate  the  beam  shift,  the  registration  error  decreased  by  88%  in  the  a  direction  and 
45%  in  the  /?  direction. 


Table  4  Mean  error  in  the  shift  estimates 


Parameter 

Mean  Error 
(pixels) 

Mean  Error 
15  th  Frame 
Removed 

Decrease  in 
Error 

a  -  beam  tracking  on 

-0.08 

-0.07 

88% 

a  -  beam  tracking  off 

-1.07 

-0.60 

(3  -  beam  tracking  on 

0.16 

0.16 

45% 

(3  -  beam  tracking  off 

0.36 

0.29 

y  -  beam  tracking  on 

-1.57 

-1.48 

b  -  beam  tracking  on 

1.36 

1.34 
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4.2  Measured  Data 


Using  the  setup  described  in  Figure  15  and  Table  2,  the  algorithm  was  used  to 
estimate  the  shifts  and  target  using  the  same  parameters  from  the  simulated  data  and 
shown  in  Table  3.  The  estimated  images  produced  by  the  algorithm  are  shown  in  Figure 
24  and  Figure  25.  Visually  comparing  Figure  24  and  Figure  25,  the  tracked  beam  case  in 
Figure  24  is  clearly  better  resolved  than  the  untracked  case  in  Figure  25.  Additionally, 
when  the  beam  is  tracked  the  estimated  image  is  slightly  larger  than  the  image  estimated 
when  the  beam  is  not  tracked.  This  is  shown  more  clearly  in  Figure  26  which  is  looking 
only  at  the  fringes  of  the  estimate  image.  There  is  clearly  more  information  on  the 
fringes  when  the  beam  is  being  tracked. 


Figure  24:  Estimated  target  using  measured  data  -  beam  tracking  on 
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Figure  25:  Estimated  target  using  measured  data  -  beam  tracking  off 


Beam  On 


Beam  Off 


Figure  26:  Difference  in  estimated  target  at  the  fringes 


The  estimated  shifts,  given  in  Appendix  B,  are  summarized  by  using  Equation  36 
to  compute  the  error  at  each  frame  when  compared  to  the  true  shift  and  is  shown  in 
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Figure  27,  Figure  28,  Figure  29  and  Figure  30.  The  true  shifts,  controlled  in  the 
MATLAB  environment  are  the  same  shifts  used  in  the  simulated  data.  However,  they 
need  to  be  scaled  by  1.104  as  was  found  in  Equation  35  to  obtain  the  true  shifts  in  the 
measured  data.  Similar  to  the  experimental  data,  frame  15  appears  to  be  an  outlier  and 
the  mean  is  also  calculated  with  that  frame  removed  and  is  also  summarized  in  Table  5. 


Figure  27:  Error  in  shift  estimate  for  a  using  measured  data 
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Shift  Error  (pixel) 


Figure  28:  Error  in  shift  estimate  for  p  using  measured  data 


Figure  29:  Error  in  shift  estimate  for  y  using  measured  data 
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Figure  30:  Error  in  shift  estimate  for  £  using  measured  data 


Table  5  Mean  error  in  the  shift  estimates  using  measured  data 


Parameter 

Mean  Error 
(pixels) 

Mean  Error 
15  th  Frame 
Removed 

Decrease  in 
Error 

a  -  beam  tracking  on 

-0.46 

-0.39 

55% 

a  -  beam  tracking  off 

-0.96 

-0.86 

[3  -  beam  tracking  on 

0.06 

0.04 

78% 

P  -  beam  tracking  off 

0.21 

0.18 

y  -  beam  tracking  on 

-1.55 

-1.45 

b  -  beam  tracking  on 

2.05 

1.98 

Analyzing  the  information  provided  in  the  previous  graphs  show  that  in  most 
frames  the  algorithm  provides  a  better  estimate  of  the  shifts  when  the  beam  is  tracked. 
Looking  at  the  mean  error  with  the  15th  frame  removed  shows  a  55%  and  78%  reduction 
in  the  a  and  /?  shift  error  respectively  when  the  beam  is  tracking.  Thus  a  lower  RMSE 
and  higher  resolution  image  would  be  expected  when  the  beam  tracking  is  on  and  the 
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images  are  averaged.  Due  to  the  scaling  from  the  camera  and  uncertainties  in  measured 
data,  a  RMSE  of  the  photon  error  cannot  be  calculated  since  the  true  image  is  not  known. 

4.3  Analysis  of  Results 

Examining  both  measured  and  simulated  data,  the  proposed  algorithm  provides  an 
improvement  in  registration  performance  when  the  beam  is  tracked.  The  reduction  in 
shift  error  is  similar  between  the  simulated  and  experimental  data  sets.  This  is  significant 
due  to  the  additional  challenges  associated  with  using  measured  data.  Specifically,  the 
PSF  and  the  summation  of  the  beam  shifts  from  Equation  31  cannot  be  known  exactly 
when  working  with  measured  data.  Mathematically,  the  beam  sum  term  can  be 
simplified  since  the  summation  of  numerous  Gaussian  beams  that  are  shifted  results  in 
another  Gaussian  beam  that  has  a  larger  standard  deviation  [6].  However  this 
mathematical  simplification  will  not  be  equal  to  the  exact  sum  of  the  shifted  beams  when 
there  are  a  limited  number  of  frames  of  data  as  is  the  case  in  these  data  sets.  These 
unknown  factors  along  with  the  noise  introduced  in  the  measured  data  create  a  situation 
that  is  not  as  ideal  as  working  in  a  purely  simulated  environment. 

A  stopping  criterion  for  the  beam  and  global  shift  estimates  and  image 
deconvolution  was  necessary  in  this  experiment.  It  was  observed  that  once  the  estimated 
shifts  were  identical  to  the  previous  iteration,  the  estimates  would  not  further  improve  but 
instead  slowly  diverge  from  the  solution.  This  was  due  to  the  algorithm’s  design  in  that 
at  each  iteration  it  attempts  to  make  the  estimated  image  look  more  like  the  detected 
image  which  includes  noise.  Thus,  a  stopping  criterion  was  set  that  once  the  estimated 
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shifts  didn’t  change  from  the  previous  estimates,  the  algorithm  stopped  updating  the  shift 
estimates. 

Of  additional  interest  in  the  results  was  the  ability  of  the  algorithm  to  create  a 
larger  image  when  the  beam  is  tracked.  This  resulted  in  more  information  on  the  target 
being  estimated.  This  is  a  direct  result  of  tracking  the  beam  at  each  frame  instead  of 
considering  that  the  beam  is  stationary.  As  the  amount  of  beam  wander  increases,  this 
advantage  proportionally  increases  however  the  algorithm’s  error  will  increase  if  there  is 
too  much  beam  wander  to  accurately  register  frames. 

To  further  validate  the  algorithm’s  increased  performance  abilities,  the  shift 
estimates  from  the  proposed  algorithm  are  compared  to  the  estimates  obtained  from  a 
cross-correlation  registration  algorithm  [2],  The  results  are  shown  in  Figure  31,  Figure 
32,  Figure  33  and  Figure  34  show  the  significant  improvement  in  registration 
performance  the  proposed  algorithm  provides  over  a  cross -correlation  algorithm  used  for 
image  registration  using  both  simulated  and  measured  data.  Similar  to  the  previous 
results,  frame  15  could  be  considered  an  outlier  and  thus  the  mean  error  with  that  frame 
removed  is  shown  in  Table  6  using  the  simulated  data  and  in  Table  7  using  the 
experimental  data.  Both  data  sets  show  a  19%  to  83%  reduction  in  shift  registration  error 
when  using  the  proposed  algorithm  to  register  frames. 
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Figure  31:  Error  in  shift  estimate  for  a  using  simulated  data 


Figure  32:  Error  in  shift  estimate  for  p  using  simulated  data 
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Figure  33:  Error  in  shift  estimate  for  a  using  measured  data 


Figure  34:  Error  in  shift  estimate  for  p  using  measured  data 
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Table  6  Mean  error  in  the  shift  estimates  using  simulated  data 


Parameter 

Mean  Error 
(pixels) 

15  th  Frame 
Removed 
Mean  Error 
(pixels) 

Decrease  in 
Error 

a  -  proposed  algorithm 

-0.14 

-0.11 

79% 

a  -  cross-correlation 

-1.82 

-0.52 

(3  -  proposed  algorithm 

-0.11 

-0.09 

43% 

(3  -  cross-correlation 

0.55 

0.16 

Table  7  Mean  error  in  the  shift  estimates  using  measured  data 


Parameter 

Mean  Error 
(pixels) 

15  th  Frame 
Removed 
Mean  Error 
(pixels) 

Decrease  in 
Error 

a  -  proposed  algorithm 

-0.46 

-0.39 

19% 

a  -  cross-correlation 

-1.76 

-0.48 

[3  -  proposed  algorithm 

0.06 

-0.04 

83% 

P  -  cross-correlation 

0.61 

0.23 
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V.  Conclusions  and  Recommendations 


This  section  details  the  conclusions  that  can  be  made  from  the  results  and  analysis 
of  this  research.  Additionally,  future  related  follow-on  research  to  this  effort  is  presented. 

5.1  Conclusions  of  Research 

The  results  from  this  research  prove  that  under  certain  circumstances  beam 
wander  caused  by  atmospheric  turbulence  can  be  tracked  independently  of  scene  shifting. 
This  results  in  a  superior  registered  image  after  post-processing  as  well  as  tracking 
information  for  the  beam  and  scene  separately.  The  algorithm  proposed  in  this  research 
was  proven  using  both  simulated  and  a  measured  data  set  to  provide  an  improved 
performance  when  compared  to  a  cross-correlation  standard  image  registration  algorithm 
which  does  not  track  beam  wander  but  considers  it  to  be  stationary. 

Unique  to  this  research  was  the  use  of  a  hybrid  approach  to  collecting  measured 
data.  Data  was  displayed  on  a  computer  screen  so  that  the  true  shifts  in  the  beam  and 
scene  could  be  controlled  and  known.  A  camera  was  used  to  capture  the  scene  displayed 
on  the  screen  thus  unknown  parameters  associated  with  true  measured  data  such  as  noise 
are  intact.  This  hybrid  approach  allowed  measured  data  to  be  captured  in  a  space  limited 
laboratory  environment  that  would  have  taken  several  kilometers  on  a  test  range  to 
collect. 

The  capabilities  of  this  algorithm  have  potential  significant  defense  applications. 
The  ability  to  reduce  registration  error  results  in  an  image  that  will  have  greater 
resolution  providing  the  end  user  with  more  information  on  the  target.  The  defense  and 
intelligence  applications  could  include  damage  assessment  of  a  laser  weapon  strike  or 
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more  detailed  information  on  a  target  in  space  situational  awareness  and  LADAR 
imaging  scenarios.  Specifically,  the  ability  of  the  algorithm  to  provide  greater  scene 
tracking  capability  in  the  presence  of  beam  wander  could  have  applications  to  improving 
the  performance  of  the  USAF’s  airborne  laser  weapon  system. 

5.2  Recommendations  for  Future  Research 

The  algorithm  provides  an  improved  performance  in  registration  but  just  how 
significant  is  this  improvement  to  various  applications  is  not  known.  An  analysis  on  the 
capabilities  this  improved  performance  could  provide  is  a  topic  that  could  be  addressed. 
Additionally,  this  algorithm  as  written  does  not  provide  real  time  feedback.  In  future 
work,  the  processing  time  could  possibly  be  improved  to  provide  near  real  time  tracking 
information.  This  could  result  in  significantly  improving  beam  and  scene  tracking 
information  for  defense  and  commercial  applications. 

Future  follow-on  work  to  this  research  could  emphasize  reducing  the  restraints  on 
the  assumptions  made  to  scope  the  level  of  work.  Currently,  the  algorithm  only  has  the 
ability  to  work  with  a  fixed  and  known  PSF.  A  blind  deconvolution  approach  to  estimate 
the  PSF  at  each  iteration  would  improve  the  performance  and  applications  of  the 
algorithm.  Another  factor  is  that  different  types  of  registration  error  are  not  included  in 
this  algorithm  and  are  both  possible  and  likely.  These  include  scaling  or  off  axis  rotation 
of  the  scene  between  images.  Future  work  that  reduces  or  eliminates  these  constraints 
would  improve  performance  and  broaden  its  relevance. 

Lastly,  future  work  could  be  done  on  using  true  measured  data  from  a  test  range 
and  testing  the  algorithm’s  ability  to  prove  it  can  provide  superior  performance  with  this 
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data.  The  hybrid  approach  to  collecting  the  measured  data  used  in  this  research  may  find 
skeptics  who  are  leery  that  the  algorithm  would  not  perform  as  well  if  true  measured  data 
was  to  be  used. 
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Appendix  A 


This  appendix  provides  a  detailed  derivation  of  the  first  term  in  the  conditional 
expectation  of  the  log-likelihood  function  in  Chapter  111. 

A.l  Conditional  Expectation  of  the  Complete  Data  Log  Likelihood 

The  expression  shown  below  in  Equation  38  was  the  original  equation  to  be 
simplified  and  solved  for  in  deriving  of  the  conditional  expectation  of  the  log-likelihood 
function  in  Chapter  111. 

E[5k(x,z)  ln(o(z)ft(z  —  yf)h(x  —  z  —  ak))  |  oold,  dfc(x),<w,KfcW]  (38) 

First,  two  statistically  independent  Poisson  random  variables,  dt and  d2,  are 
defined  as  shown  in  Equation  39  and  Equation  40.  The  new  Poisson  random  variable,  dlt 
is  one  frame  of  the  complete  data  at  a  single  point  z0.  The  Poisson  random  variable,  d2, 
is  the  sum  of  all  the  frames  and  pixel  points  in  the  complete  data  except  for  the  point 
in  d1?  the  Poisson  background  noise  is  also  added  to  this  sum. 

d1  =  d(x,zQ)  (39) 

K  N 

d{x,z)  —  d(x,  z0)  +  B(x)  (4Q) 

k  z 

These  two  random  variables,  djand  d2,  are  defined  to  have  some  mean,  m,  as 
shown  in  Equation  41  and  Equation  42. 


EldjJ  =  m1 

(41) 

E  [d2]  =  m2 

(42) 
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Looking  back  at  Equation  10,  the  incomplete  data  is  related  to  dx and  d2  as  shown 


in  Equation  43. 


d(x)  = 


K  N 

cLk(x,z )  +  B(x )  =  d1  +  d2 


(43) 


Using  the  Poisson  PMF,  the  joint  probability  of  the  random  variables  dxand  d2 
with  means  m^nd  m2  is  given  in  Equation  44. 


P(di,  d2) 


(44) 


The  variable  d2  can  be  removed  from  Equation  44  to  allow  simplification  using 
the  equality  in  Equation  45.  The  new  joint  probability  is  shown  in  Equation  46. 


d2  —  d  d^ 


(45) 


P(d,  dU 


/mdle~mi\  fmd2-dle-mt 

l  dd  )  V  (d-dj! 


(46) 


Using  Bayes  theorem  [6],  the  conditional  expectation  of  dx  is  found  and 


simplified  in  Equation  47. 


Pldjd) 


P(d,  dx) 

P  (d) 


fmd  dle  m2\ 

{  (d-dj!  J 


(m1  +  m2)de~(m  i+m2) 
d! 


(dljm^1  md  dl 
(irij  +  m2)d(d1\)  (d  —  djJ! 


(dljm^1  dl 

(tti-l  +  77i2)(di+d_£ii)(d1!)  (d  -  dx) ! 

(d-L !)  (d  —  d^! 


(47) 


The  final  form  of  the  conditional  expectation  in  Equation  47  is  similar  to  the  PMF 
of  a  binomial  random  variable  shown  in  Equation  48.  The  binomial  PMF  describes  the 
probability  of  getting  exactly  n  successes  in  k  trials  for  an  event  with  a  probability  of 
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success  p  [6].  The  mean  of  a  binomial  is  shown  in  Equation  49.  The  relationships 


between  Equation  47  and  the  binomial  PMF  are  shown  in  Table  A.l. 


P(n,  p)  = 


n! 


k\  (n  —  fc)! 


pfc(l  -p) 


n-k 


(48) 


E[Binominal(n,p)]  =  np 


(49) 


Table  A.l  Relationship  between  the  conditional  expectation  log  likelihood  function  and  the  binomial 

PMF 


Number  of  trials 

d 

Number  of  success 

dj 

Probability  of  success 

m1 

(m1  +  m2) 

Probability  of  failure 

m2 

( m1  +m2) 

Mean 

d(  mi  ) 
\m1  +  m2J 

Converting  back  to  the  original  notation  in  Equation  38,  the  conditional 
expectation  of  the  complete  data  log-likelihood  is  shown  in  Equation  50. 


E  [dk{x,z)  ln(o(»i>0  -  yk)h(x  -z-  ak ))  |  6oW,dk(x),<w,KfcW] 

=  E[4(x,z)|ooW,  dk(x),  a£ld,  ykld]  In (o(z)b(z  -  yk)l i(x  -  z  -  ak)) 


cL(x)dold (z)b{z  -  ykld)h(x  -  z  -  akld ) 
ikol\x) 


In  (o(z)i)(z  —  —  z  —  afc)) 


Where  shown  again  for  easy  reference, 

N 

l(x)  —  ^  dold(z,w)b(z  —  y£ld)h(x  —  z  —  akld ). 


old  ( 


(50) 


(51) 
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Appendix  B 


This  appendix  provides  the  true  shifts  in  the  MATLAB  environment,  as  well  as 
the  shifts  estimated  by  the  algorithm  using  both  the  simulated  data  and  measured  data. 

B.l  Simulated  Data  Results 

The  true  shift  and  the  estimated  shifts  from  the  algorithm  when  using  simulated 
data  are  given  explicitly  in  Table  B.l,  B.2,  B.3  and  B.4. 

B.2  Measured  Data  Results 

The  scaled  true  shifts  and  the  estimated  shifts  from  the  algorithm  when  using 
measured  data  are  given  explicitly  in  Table  B.5,  B.6,  B.7  and  B.8. 
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